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Abstract 


Our paper computationally explores Ricker’s predator satiation model with the objective 
of studying how the extinction dynamics of an animal species having a two-stage life- 
cycle is affected by a sudden spike in mortality due to an extraneous extreme event. Our 
simulation model has been designed and implemented using sockeye salmon population 
data based on a stochastic version of Ricker’s model; with the shock size being reflected 
by a sudden reduction in the carrying capacity of the environment for this species. Our 
results show that even for a relatively marginal increase in the negative impact of an 
extreme event on the carrying capacity of the environment, a species with an otherwise 
stable population may be driven close to extinction. 


Key words: Ricker’s model, extinction dynamics, extreme event, Monte Carlo 
simulation 


Background and research objective 


PVA approaches do not normally consider the risk of catastrophic extreme events under 
the pretext that no population size can be large enough to guarantee survival of a species 
in the event of a large-scale natural catastrophe. ''! Nevertheless, it is only very intuitive 
that some species are more “delicate” than others; and although not presently under any 
clearly observed threat, could become threatened with extinction very quickly if an 
extreme event was to occur even on a low-to-moderate scale. The term “extreme event” is 
preferred to “catastrophe” because catastrophe usually implies a natural event whereas; 
quite clearly; the chance of man-caused extreme events poses a much greater threat at 
present to a number of animal species as compared to any large-scale natural catastrophe. 
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An animal has a two-stage life cycle when; in the first stage, newborns become immature 
youths and in the second stage; the immature youths become mature adults. Therefore, in 
terms of the stage-specific approach, if Y; denotes the number of immature young in 
stage t and A; denotes the number of mature adults, then the number of adults in year t + 
1 will be some proportion of the young, specifically those that survive to the next 
(reproductive) stage. Then the formal relationship between the number of mature adults 
in the next stage and the number of immature youths at present may be written as 
follows: 


At+1=AYt 


Here A is the survival probability, i.e. it is the probability of survival of a youth to 
maturity. The number of young next year will depend on the number of adults in t: 


Yi+1 =f (Ad) 
Here f describes the reproduction relation between mature adults and next year’s young. 


This is a straightforward system of simultaneous difference equations which may be 
analytically solved using a variation of the cobwebbing approach. ®! The solution process 
begins with an initial point (Y;, Ai) and iteratively determines the next point (Y2, Az). If 
predator satiation is built into the process, then we simply end up with Ricker’s model: 


—At/K 
Yi+1=aAce 


Here a is the maximum reproduction rate (for an initial small population) and K is the 
population size at which the reproduction rate is approximately half its maximum "! 
Putting £ = '/x we can re-write Ricker’s equation as follows: 


Yt+1 =aAe He 


It has been shown that if (Yo, Ao) lies within the first of three possible ranges, (Y;, An) 
approaches (0, 0) in successive years and the population becomes extinct. If (Yo, Ao) lies 
within the third range then (Yn, An) equilibrate to a steady-state value of (Y*, A*). 
Populations that begin with (Yo, Ao) within the second range oscillate between (Y*, 0) 
and (0, A*). Such alternating behavior indicates one of the year classes, or cohorts, 
become extinct while the other persists i.e. adult breeding stock appear only every other 
year. Thus the model reveals that three quite different results occur depending initially 
only on the starting sizes of the population and its distribution among the two stages. "“! 


We use the same basic model in our research but instead of analytically solving the 
system of difference equations, we use the same to simulate the population dynamics as a 
stochastic process implemented on an MS-Excel spreadsheet. Rather than using a closed- 
form equation like Ricker’s model to represent the functional relationship between Y; + | 
and A;, we use a Monte Carlo method to simulate the stage-transition process within 
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Ricker’s framework; introducing a massive perturbation with a very small probability in 
order to emulate a catastrophic event.! 


Conceptual framework 


We have a formulated a stochastic population growth model with an inbuilt capacity to 
generate an extreme event based on a theoretical probability distribution. The non- 
stochastic part of the model corresponds to Ricker’s relationship between Y;, +, and Aj. 
The stochastic part has to do with whether or not an extreme event occurs at a particular 
time point. The gamma distribution has been chosen to make the probability distribution 
for the extreme event a skewed one as it is likely to be in reality. Instead of analytically 
solving the system of simultaneous difference equations iteratively in some variation of 
the cobwebbing method, we have used them in a spreadsheet model to simulate the 
population growth over a span of ten time periods. 


We apply a computational methodology whereby the initial number of immature young is 
hypothesized to either attain the expected number predicted by Ricker’s model or 
drastically fall below that number at the end of every stage, depending on whether an 
extraneous extreme event does not occur or actually occurs. The mortalities as a result of 
an extreme event at any time point is expressed as a percentage of the pristine population 
size for a clearer comparative view. 


Model building 


Among various faunal species, the population dynamics of the sockeye salmon 
(oncorhynchus nerka) has been most extensively studied using Rickert’s model. Salmon 
are unique in that they breed in particular fresh water systems before they die. Their 
offspring migrates to the ocean and upon reproductive maturity, they are guided by a 
hitherto unaccounted instinctive drive to swim back to the very same fresh waters where 
they were born to spawn their own offspring and perish. Salmon populations thus are 
very sensitive to habitat changes and human activities that have a negative impact on 
riparian ecosystems that serve as breeding grounds for salmon can adversely affect the 
peculiar life-cycle of the salmon. Many of the ancient salmon runs (notably those in 
California river systems) have now gone extinct and it is our hypothesis that an even 
seemingly stabilized population can be rapidly driven to extinction due to the effect of an 
extraneous (quite possibly man-made) extreme event with the capacity to cause mass 
mortality. The following table shows the four-year averages of the sockeye salmon 
population in the Skeena river system in British Columbia in the first half of the twentieth 
century. 
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Year | Population (in thousands) 
1908 1,098 
1912 740 
1916 714 
1920 615 
1924 706 
1928 510 
1932 278 
1936 448 
1940 528 
1944 639 
1948 523 


(Source: http://www-rohan.sdsu.edu/~jmahaffv/courses/s00/math 12 INectures/product_rule/product.html#Ricker'sModel) 


A non-linear least squares best-fit to Ricker’s model is obtained for the above set of data 
is obtained as follows: 


Minimize e° = \'[d, — {ad e+]? , where d, is the actual population size in year t. 
t t pop y 


t=1 


The necessary conditions to the above least squares best-fit problem is obtained as 
follows: 


A(e’)/Oa = O(e”)/OB = 0; whereby we get a* ~ 1.54 and B* ~7.8 x 10° 


Plugging these parameters into Ricker’s model indeed yields a fairly good approximation 
of the salmon population stabilization in the Skeena river system in the first half of the 
previous century. 


As the probability distribution of an extraneous extreme event is likely to be a highly 
skewed one, we have generated our random variables from the cumulative distribution 
function (cdf) of the gamma distribution rather than the normal distribution. The 
distribution boundaries are fixed by generating random integers in the range 1 to 100 and 
using these random integers to define the shape and scale parameters of the gamma 
distribution. The gamma distribution performs better than the normal distribution when 
the distribution to be matched is highly right-skewed; as is desired in our model. The 
combination of a large variance and a lower limit at zero makes fitting a normal 
distribution rather unsuitable in such cases.'°! The probability density function of the 
gamma distribution is given as follows: 


f(x, a,b) = {b°T(a)}'x“'e*” forx>0 
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Here a > 0 is the shape parameter and B > 0 is the scale parameter of the gamma 
distribution. The cumulative distribution function may be expressed in terms of the 
incomplete gamma function as follows: 


F (x, a, b) = [fend =y(a,x/b)/T(a) 


In our spreadsheet model, we have F (R, R/2, 2) as our cdf of the gamma distribution. 
Here R is an integer randomly sampled from the range 1 to 100. An interesting statistical 
result of having these values for x, a and f is that the cumulative gamma distribution 
value becomes equalized with the value [1 - y? (R)] having R degrees of freedom, thus 
allowing y” goodness-of-fit tests. ! 


Our model is specifically designed to simulate the extinction dynamics of sockeye 
salmon population using a stochastic version of Ricker’s model; with the shock size 
being based on a sudden reduction in the parameter K i.e. the carrying capacity of the 
environment for this species. The model parameters are same as those of Ricker’s model 
i.e. a and B (which is the reciprocal of K). We have kept a constant at all times at 1.54, 
which was the least squares best-fit value obtained for that parameter. We have kept a B 
of 0.00078 (i.e. the best-fit value) when no extreme event occurs and have varied the B 
between 0.00195 and 0.0156 (i.e. between 2.5 times to 20 times the best-fit value) for 
cases where an extreme event occurred. We have a third parameter c which is basically a 
‘switching constant’ that determines whether an extreme event occurs or not. The switch 
is turned on triggering an extreme event when a random draw from a cumulative gamma 
distribution yields a value less than or equal to c. Using F (R, R/2, 2) as our cdf of the 
gamma distribution where R is a randomly drawn integer in the range (1, 100) means that 
the cumulative gamma function will randomly select from the approximate interval 0.518 
~ 0.683. By fixing the value of c at 0.5189 in our model we have effectively reduced the 
probability of occurrence of an extreme event to a miniscule magnitude relative to that of 
an extreme event not occurring. We have used the sockeye salmon population data from 
the table presented earlier For each level of the B parameter, we simulated the system and 
observed the maximum possible number of mortalities from an extreme event at that 
level of B. The results are reported below. 
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Results obtained from the simulation model 


Worst-case effect of extreme event on sockeye salmon population 
20% co 
15% | 
10% | 


5% 
O% 


population) 


Surviving population size 
(in terms of % of pristine 


0 0.005 0.01 0.015 0.02 


Shock size (in terms of impact on carrying capacity) 


We made 100 independent simulation runs for each of the eight levels of 6. The low 
probability of extreme event assigned in our study yielded a mean of 1.375 for the 
number of observed worst-case scenarios (i.e. situations of maximum mortality) with a 
standard deviation of approximately 0.92. The worst-case scenarios for our choice of 
parameters necessarily occur if the extreme event occurs in the first time point when the 
species population is at its maximum size. Our model shows that in worst-case scenarios, 
the size of surviving population after an extreme event that could seed the ultimate 
recovery of the species to pre-catastrophe numbers (staying within the broad framework 
of Ricker’s model) drops from about 18% of the pristine population size for a shock size 
corresponding to 2.5 times the best-fit B; to only about 0.000005% of the pristine 
population size for a shock size corresponding to 20 times the best-fit p. 


Therefore, if the minimum required size of the surviving population is at least say 20% of 
the pristine population in order to survive and recover to pre-catastrophe numbers, the 
species could go extinct if an extreme event caused a little more than two-fold decrease in 
the environmental carrying capacity! Even if the minimum required size for recovery was 
relatively low at say around 2% of the pristine population, an extreme event that caused a 
five-fold decrease in the environmental carrying capacity could very easily force the 
species to the brink of extinction. An immediate course of future extension of our work 
would be allowing the fecundity parameter a to be affected by extreme events as is very 
likely in case of say a large-scale chemical contamination of an ecosystem due to a faulty 
industrial waste-treatment facility. 
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Conclusion 


Our study has shown that even for a relatively marginal 2.5-fold decrease in the 
environmental carrying capacity due to an extreme event, a worst-case scenario could 
mean a mortality figure well above 80% of the pristine population. As a guide for future 
PVA studies we may suggest that one should not be deterred simply by the notion that 
extreme events are uncontrollable and hence outside the purview of computational 
modeling. Indeed the effect of an extreme event can almost always prove to be fatal for a 
species but nevertheless, as our study shows, there is ample scope and justification for 
future scientific enquiries into the relationship between survival probability of a species 
and the adverse impact of an extreme event on ecological sustainability. 
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